/* Bessel J_0(x) function in double precision */

#include <math.h>

double dbesj0(double x)
{
    int k;
    double w, t, y, v, theta;
    static double a[8] = {
        -2.3655394e-12, 4.70889868e-10, 
        -6.78167892231e-8, 6.7816840038636e-6, 
        -4.340277777716935e-4, 0.0156249999999992397, 
        -0.2499999999999999638, 0.9999999999999999997
    };
    static double b[65] = {
        6.26681117e-11, -2.2270614428e-9, 
        6.62981656302e-8, -1.6268486502196e-6, 
        3.21978384111685e-5, -5.00523773331583e-4, 
        0.0059060313537449816, -0.0505265323740109701, 
        0.2936432097610503985, -1.0482565081091638637, 
        1.9181123286040428113, -1.13191994752217001, 
        -0.1965480952704682, 
        4.57457332e-11, -1.5814772025e-9, 
        4.55487446311e-8, -1.0735201286233e-6, 
        2.02015179970014e-5, -2.942392368203808e-4, 
        0.0031801987726150648, -0.0239875209742846362, 
        0.1141447698973777641, -0.2766726722823530233, 
        0.1088620480970941648, 0.5136514645381999197, 
        -0.2100594022073706033, 
        3.31366618e-11, -1.1119090229e-9, 
        3.08823040363e-8, -6.956602653104e-7, 
        1.23499947481762e-5, -1.66295194539618e-4, 
        0.0016048663165678412, -0.0100785479932760966, 
        0.0328996815223415274, -0.0056168761733860688, 
        -0.2341096400274429386, 0.2551729256776404262, 
        0.2288438186148935667, 
        2.38007203e-11, -7.731046439e-10, 
        2.06237001152e-8, -4.412291442285e-7, 
        7.3107766249655e-6, -8.91749801028666e-5, 
        7.34165451384135e-4, -0.0033303085445352071, 
        0.0015425853045205717, 0.0521100583113136379, 
        -0.1334447768979217815, -0.1401330292364750968, 
        0.2685616168804818919, 
        1.6935595e-11, -5.308092192e-10, 
        1.35323005576e-8, -2.726650587978e-7, 
        4.151324014176e-6, -4.43353052220157e-5, 
        2.815740758993879e-4, -4.393235121629007e-4, 
        -0.0067573531105799347, 0.0369141914660130814, 
        0.0081673361942996237, -0.257338128589888186, 
        0.0459580257102978932
    };
    static double c[70] = {
        -3.009451757e-11, -1.4958003844e-10, 
        5.06854544776e-9, 1.863564222012e-8, 
        -6.0304249068078e-7, -1.47686259937403e-6, 
        4.714331342682714e-5, 6.286305481740818e-5, 
        -0.00214137170594124344, -8.9157336676889788e-4, 
        0.04508258728666024989, -0.00490362805828762224, 
        -0.27312196367405374426, 0.04193925184293450356, 
        -7.1245356e-12, -4.1170814825e-10, 
        1.38012624364e-9, 5.704447670683e-8, 
        -1.9026363528842e-7, -5.33925032409729e-6, 
        1.736064885538091e-5, 3.0692619152608375e-4, 
        -9.2598938200644367e-4, -0.00917934265960017663, 
        0.02287952522866389076, 0.10545197546252853195, 
        -0.16126443075752985095, -0.19392874768742235538, 
        2.128344556e-11, -3.1053910272e-10, 
        -3.34979293158e-9, 4.50723289505e-8, 
        3.6437959146427e-7, -4.46421436266678e-6, 
        -2.523429344576552e-5, 2.7519882931758163e-4, 
        9.7185076358599358e-4, -0.00898326746345390692, 
        -0.01665959196063987584, 0.11456933464891967814, 
        0.07885001422733148815, -0.23664819446234712621, 
        3.035295055e-11, 5.486066835e-11, 
        -5.01026824811e-9, -5.0124684786e-9, 
        5.8012340163034e-7, 1.6788922416169e-7, 
        -4.373270270147275e-5, 1.183898532719802e-5, 
        0.00189863342862291449, -0.0011375924956163613, 
        -0.03846797195329871681, 0.02389746880951420335, 
        0.22837862066532347461, -0.06765394811166522844, 
        1.279875977e-11, 3.5925958103e-10, 
        -2.28037105967e-9, -4.852770517176e-8, 
        2.8696428000189e-7, 4.40131125178642e-6, 
        -2.366617753349105e-5, -2.4412456252884129e-4, 
        0.00113028178539430542, 0.0070847051391978908, 
        -0.02526914792327618386, -0.08006137953480093426, 
        0.16548380461475971846, 0.14688405470042110229
    };
    static double d[52] = {
        1.059601355592185731e-14, -2.71150591218550377e-13, 
        8.6514809056201638e-12, -4.6264028554286627e-10, 
        5.0815403835647104e-8, -1.76722552048141208e-5, 
        0.16286750396763997378, 2.949651820598278873e-13, 
        -8.818215611676125741e-12, 3.571119876162253451e-10, 
        -2.63192412099371706e-8, 4.709502795656698909e-6, 
        -0.005208333333333283282, 
        7.18344107717531977e-15, -2.51623725588410308e-13, 
        8.6017784918920604e-12, -4.6256876614290359e-10, 
        5.0815343220437937e-8, -1.7672255176494197e-5, 
        0.16286750396763433767, 2.2327570859680094777e-13, 
        -8.464594853517051292e-12, 3.563766464349055183e-10, 
        -2.631843986737892965e-8, 4.70950234228865941e-6, 
        -0.0052083333332278466225, 
        5.15413392842889366e-15, -2.27740238380640162e-13, 
        8.4827767197609014e-12, -4.6224753682737618e-10, 
        5.0814848128929134e-8, -1.7672254763876748e-5, 
        0.16286750396748926663, 1.7316195320192170887e-13, 
        -7.971122772293919646e-12, 3.544039469911895749e-10, 
        -2.631443902081701081e-8, 4.709498228695400603e-6, 
        -0.005208333331514365361, 
        3.84653681453798517e-15, -2.04464520778789011e-13, 
        8.3089298605177838e-12, -4.6155016158412096e-10, 
        5.081326369646665e-8, -1.76722528311426167e-5, 
        0.1628675039665006593, 1.3797879972460878797e-13, 
        -7.448089381011684812e-12, 3.51273379710695978e-10, 
        -2.630500895563592722e-8, 4.709483934775839193e-6, 
        -0.0052083333227940760113
    };

    w = fabs(x);
    if (w < 1) {
        t = w * w;
        y = ((((((a[0] * t + a[1]) * t + 
            a[2]) * t + a[3]) * t + a[4]) * t + 
            a[5]) * t + a[6]) * t + a[7];
    } else if (w < 8.5) {
        t = w * w * 0.0625;
        k = (int) t;
        t -= k + 0.5;
        k *= 13;
        y = (((((((((((b[k] * t + b[k + 1]) * t + 
            b[k + 2]) * t + b[k + 3]) * t + b[k + 4]) * t + 
            b[k + 5]) * t + b[k + 6]) * t + b[k + 7]) * t + 
            b[k + 8]) * t + b[k + 9]) * t + b[k + 10]) * t + 
            b[k + 11]) * t + b[k + 12];
    } else if (w < 12.5) {
        k = (int) w;
        t = w - (k + 0.5);
        k = 14 * (k - 8);
        y = ((((((((((((c[k] * t + c[k + 1]) * t + 
            c[k + 2]) * t + c[k + 3]) * t + c[k + 4]) * t + 
            c[k + 5]) * t + c[k + 6]) * t + c[k + 7]) * t + 
            c[k + 8]) * t + c[k + 9]) * t + c[k + 10]) * t + 
            c[k + 11]) * t + c[k + 12]) * t + c[k + 13];
    } else {
        v = 24 / w;
        t = v * v;
        k = 13 * ((int) t);
        y = ((((((d[k] * t + d[k + 1]) * t + 
            d[k + 2]) * t + d[k + 3]) * t + d[k + 4]) * t + 
            d[k + 5]) * t + d[k + 6]) * sqrt(v);
        theta = (((((d[k + 7] * t + d[k + 8]) * t + 
            d[k + 9]) * t + d[k + 10]) * t + d[k + 11]) * t + 
            d[k + 12]) * v - 0.78539816339744830962;
        y *= cos(w + theta);
    }
    return y;
}

